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Heating and Cooling Protostellar Disks 

S. Hirose 1 and N. J. Turner 2 

ABSTRACT 

We examine heating and cooling in protostellar disks using 3-D radiation- 
MHD calculations of a patch of the Solar nebula at 1 AU, employing the shearing- 
box and flux-limited radiation diffusion approximations. The disk atmosphere is 
ionized by stellar X-rays, well-coupled to magnetic fields, and sustains a turbu- 
lent accretion flow driven by magneto-rotational instability, while the interior is 
resistive and magnetically dead. 

The turbulent layers heat by absorbing the light from the central star and by 
dissipating the magnetic fields. They are optically-thin to their own radiation and 
cool inefficiently The optically-thick interior in contrast is heated only weakly, by 
re-emission from the atmosphere. The interior is colder than a classical viscous 
model, and isothermal. 

The magnetic fields support an extended atmosphere that absorbs the 
starlight 1.5 times higher than the hydrostatic viscous model. The disk thick- 
ness thus measures not the internal temperature, but the magnetic field strength. 
Fluctuations in the fields move the starlight-absorbing surface up and down. The 
height ranges between 13% and 24% of the radius over timescales of several orbits, 
with implications for infrared variability. 

The fields are buoyant, so the accretion heating occurs higher in the atmo- 
sphere than the stresses. The heating is localized around current sheets, caused 
by magneto-rotational instability at lower elevations and by Parker instability at 
higher elevations. Gas in the sheets is heated above the stellar irradiation temper- 
ature, even though accretion is much less than irradiation power when volume- 
averaged. The hot optically-thin current sheets might be detectable through their 
line emission. 
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Subject headings: protoplanetary disks — magnetohydrodynamics — turbulence 
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INTRODUCTION 



Measurements of the temperatures, densities and flow velocities inside protostellar 
disks are vitally important for understanding planet formation. New probes of the con- 
ditions are made possible by the detections of infrared line emission from molecules in- 
cluding H 2 0, OH, CO, COp, CoH? and HCN, in the central 3 AU of the disks around 



T Tauri stars (INaiita et all 120031 : ICarr et all 12004 ICarr fc Najital 120081 : ISalvk et all 12008 



Pontoppidan et al. 201ol ). However interpreting the line emission requires some knowledge 
of the power source. The line excitation temperatures range from a few hundred to a few 
thousand Kelvin, and can be reached in gas near 1 AU throu gh accretion heating, if the 
stresses in the disk atmosphere approach the local gas pressure (jGlassgold et al.l 120041 1 . 



The accretion stress es in T Tauri disks likely come from magneto-rotational turbulence 
( IBalbus &: Hawleylll998l ). Magnetic fields couple only to gas that is sufficiently ionized, and 
while much of the disk is too cold for ther mal ionization, the top and bottom layers ar e 
ionized by the X-rays from the young star (jlgea fc Glassgoldlll999l : lllgner fc Nelson! 120061 1 . 
The turbulent, m agnetically- active surfa ce layers are optically-thin, and thus produce line 
emission directly ( jBai fc Goodman! 120091 ) . 



Another important heat source is the light from the central star. The starlight en- 
ters the disk at a grazing angle, and so is absorbed high in the atmosphere, well above 
the surface of optical depth unity for the disk's own longer- wavelength continuum emission 
( IChiang &: Goldreichl 119971 ). The stellar illumination dominates the atmospheric heating at 
1 AU in traditional viscous models, w hile the accretion domin ates in the interior, at typical 
accretion rates around 10 -8 M yr _1 (jPullemond et al.l 120071 ) . 



Here we report first-principles radiation-MHD calculations of the heating and cooling, 
treating the conversion of gravitational to magnetic energy through magneto-rotational in- 
stability, the dissipation of the magnetic energy as heat in the gas, an d the emissi o n and 
escape of the disk radiation. The calculations build on recent work by iFlaig et al.l (120101 ) 
in including the stellar irradiation heating and allowing distinct gas and radiation tempera- 
tures. We also choose a much lower s urface density of 1 000 g cm -3 , similar to that at 1 AU 
in the minimum- mass Solar nebula (IHayashil Il98ll ). so that temperatures are too low for 
collisional ionization. The ionization is dominated by the stellar energetic photons, which 
are absorbed in the atmosphere, leaving a magnetically-decoupled interior dead zone. 
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METHODS 



We solve the equations of rad iation-MHD using t he ZEUS code ( IStone fc Normanlll992l ). 
conserving total energy following iHirose et al.l (120061 ) . The direct starlight is treated by ex- 
plicitly solving the transfer equation, while the dis k's own emission is treated using the flux- 
limited diffusion module of Turner fc Stond ( 1200 ll ). The gas and radiation energy equations 
are adva nced through tim e in fully-implicit fashion using a generalized Newton-Raphson 
method ( IStone et al.lll992l ) with multigrid solver. Radiation emission, absorption and diffu- 
sion are computed simultaneously, as needed for accuracy when two of them almost cancel. 
We use opacities for homogeneous spherical silicate grains of normal iron content, from 
Semenov et al.l (120031 ) . Opacities for the direct starlight are averaged over the stellar spec- 
trum, while for the disk's radiation we use the Planck and Rosseland means at the local 
temperature. 

The star has one-half the mass and twice the radius of our Sun, and the spectrum of 
a 4 000 K blackbody. The starlight enters the domain top and bottom at a grazing angle, 
0.05 radians ~ 3° from the horizontal. The heating is computed each timestep by solving the 
transfer equation in the horizontally-averaged density profile. We neglect the small amount 
of absorption that would occur outside the MHD domain. 

The Ohmic resistivity varies in space and time, and is found by interpolating in a pre- 
computed lookup table whose three axes are the density, te mperature and l o cal io nization 
rate. The table holds equilibrium solutions of the simplified lllgner fc Nelson! ( 120061 ) recom- 
bination network including reactions on grains (their Model 4). The grain size of 0.1 /zm is 
chosen to give almost the same geometric cross-section per unit mass as the particle size dis- 
tribution used for the opacities. The two differ by less than 15%. The ionization rate due to 
the stellar X-rays and long-lived radionuclides is c omputed from the hor izontally-averaged 
overlying mass column using the prescriptions in Turner fc Sand (120081 ). taking an X-ray 



luminosity 2 x 10 30 erg s 1 and a long-lived radionuclide ionization rate 1.4 x 10 22 s 1 at 
nominal dust abu ndance. The magneti c diffusivity is inversely proportional to the resulting 



electron fraction (IBlaes fc Balbuslll994j ) and is very high in the weakly-ionized disk interior, 



reaching 6 x 10 cm s (or if short-lived Al were present, 2 x 10 ). To avoid short 
timesteps, we cap the diffusivity, limiting it to no more than 5 x 10 16 cm 2 s _1 . This is 
below the threshold of 5 x 10 2 for diffusion to ove rcome the shearing of radial fields under 
differential rotation (eq. 1 of Turner fc Sanoll2008l ). so weak magnetic activity is expected 
in the dead zone. However this is probably unimportant for the active layers, as we obtain 
similar MHD results with a cap ten times higher. 

T Tauri disk atmospheres appear dust-depleted by f actors of IP 1 " 4 compared with the 
interstellar medium, based on their mid-infrared colors (IFurlan et al.l 120061 . 120091 ) . Atmo- 
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spheric depletion is consistent with modeling of so me protostellar disk s over wider wave - 
length ranges, incorporating spatially- resolved data (IDuchene et al.ll2003t iPinte et al.ll2008l ). 
Furthermore, the concentration of solid material in the interior is helpful for planet forma- 
tion. We therefore assume a dust-to-gas mass ratio 10 -4 and reduce all dust effects 100 times 
below nominal, including the opacities, the cross-sections for grain surface reactions, and the 
radionuclide decay ionization. The dust is well-mixed, and thermally coupled to the gas so 
they share a single temperature, which is allowed to differ from the local radiation temper- 
ature. The gas has mean molecular weight 2.3 and an ideal equation of state with 7 = 1.4. 
A density floor of 10~ 5 times the initial midplane density is applied to prevent high Alfven 
speeds from halting the calculations. We checked the results are basically unchanged with a 
density floor ten times lower. Flow velocities also are limited to no more than 10 times the 
shear speed, |Q times the domain width. 

The domain is a box of size 0.07 x 0.28 x 0.7 AU along the radial, orbital and vertical 
directions, centered in the midplane at 1 AU and divided into 32 x 64 x 32 grid cells. 



The b oundaries are shearing-periodic, periodic, and outflow, respectively. Like iFlaig et al. 



(120101 ). we assume the magnetic fields are vertical on the top and bottom boundaries. Results 
are similar using the standard ZEUS outflow boundary condition, where the fields on the 
boundaries are computed from the extrapolated E MFs. We also add ed a small artificial 
resistivity near the boundaries for numerical stability (IHirose et al.ll2009l ). The disk radiation 
on the vertical boundaries is placed in the streaming limit, with flux equal to the product of 
the radiation energy density and the speed of light. 

Initially the gas is isothermal at 280 K and in hydrostatic balance, while the magnetic 
field is the sum of an 0.02-Gauss uniform vertical component with pressure 3 x 10 5 times 
less than the midplane gas pressure, and a vertical component varying sinusoidally in the 
radial direction with maximum pressure l/270th the midplane gas pressure. Results were 
similar in a further calculation with the pressure in the sinusoidal component 10 times less. 
We al so compared an ideal-MHD version of our calculation against the results of IFlaig et al. 
(120101 ). obtaining comparable profiles of the plasma beta parameter, turbulent Mach number 
and density fluctuations, despite our lower surface density. 

Below we relate the MHD results to the classical vi scous, irradiated disk model with 
the same surface density and accretion rate. We follow iD'Alessio et al.l ( 119981 ) in solving 
the equations of vertical hydrostatic equilibrium, radiative equilibrium in the Eddington 
approximation, and thermal balance, together with a transfer equation for the illuminating 
starlight. Scattered radiation is neglected in both the viscous and MHD calculations. 
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Fig. 1. — Heating profiles in the MHD calculation (red) and the classical viscous model with 
the same surface density and accretion rate (blue). Both have contributions from stellar 
irradiation (dotted) and turbulent or viscous dissipation (solid). The turbulent dissipation 
involves both magnetic and kinetic energy. All are averaged horizontally and over time from 
50 to 200 orbits. A grey band marks the dead zone. 



3. RESULTS 



The turbulent layers on the top and bottom are optically-thin to their own radiation, 
while the magnetically-dead interior i s optically-thick. The base of the turbulen t layer, de- 
fined by unit Elsasser number v\ z /r]Q (jSano fc Inutsukall200ll ; lTurner et al.ll2007l ). has height 
0.08 AU and Rosseland mean optical depth 0.39, so the disk photosphere falls within the 
dead zone. The starlight on the other hand is absorbed in the turbulent layers, dominating 
the mean heating above 0.11 AU (Figured]). The turbulent layers are more opaque in the 
optical than the infrared, and absorb starlight better than they re-radiate, becoming several 
times hotter than the disk radiation passing through them. The horizontal averaging in the 
irradiation prescription yields mean temperatures higher than the true irradiation tempera- 
ture of 410 K found in the viscous model. Below the disk photosphere the gas and radiation 
temperatures match. The interior, in contrast, receives no direct starlight and is heated by 
the re-radiation from the hot surface layers. The cooling is almost wholly ra diative through 



out. Turbulent advection is negligible, and unlike the viscous models in iD'Alessio et al. 



(119981 ) figure 10, carries no significant heat to the interior. The heating and cooling together 
yield the temperature profile in Figure [2J 

The turbulent layers are supported by magnetic forces above 0.12 AU, while the interior 
is gas-pressure-supported. The magnetic dissipation is the second-largest contribution to the 
volume-integrated heating. The magnetic fields are buoyant and rise while dissipating, so the 
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Fig. 2. — Temperature profiles in the MHD calculation (red) and corresponding viscous 
model (blue). Shadings with the same colors indicate the penetration of the starlight to 
unit optical depth, while green shading marks the part of the disk that is opaque to its own 
radiation. The gas and dust (solid) is hotter than the disk radiation (dotted) at low optical 
depths. The profiles are averaged horizontally, and in time from 50 to 200 orbits. 



accretion heating occurs higher on average than the accretion stress (Figure |3]). By contrast, 
in the classical model the stress and heating occur at the same place and time, and both 
result from an effective viscosity, nature unspecified. The dissipation rate in our calculation 
is comparable to the local pressure p(z) times the shear rate the ratio increasing from 
0.4 near the turbulent layer base to 2 at th e boundary. Rat i o uni ty is sufficient to produce 
CO molecular line emission in the models of iGlassgold et al.l (120041 ) at lower columns around 
10 21 cm -2 . The dissipation is also high relative to the stress in two narrow layers near 
±0.07 AU where turbulent magnetic fields are destroyed on entering the dead zone. Overall, 
the stresses correspond to a mean mass flow ra te 1.3 x 10~ 8 M yr -1 , within the range among 
classical T Tauri stars (IMuzerolle et al.l 120031 ) . 

The magnetic field structure varies with height. Near the base of the turbulent layers, 
current sheets are typically close to horizontal, while high in the atmosphere, current sheets 
are most often vertical. Near the base, the current sheets separate magnetic fields moving 
radially inward and outward under the magneto-rotational instability, while higher in the at- 
mosphere, the current sheets separate magnetic fields lying along adjacent orbits and having 
out-of-phase vertical undulations. The upward slope of one field line annihilates against the 
downward slope of the next (Figure H]). The undulations come from the Parker instability, 
whose fastest-growing modes have large radial wavenumber under Keplerian differen t ial ro- 
tation. By solving in each grid cell the dispersion relation of iTerquem fc Papaloizoul (119961 ) 
eq. 33, describing localized non-axisymmetric disturbances on toroidal magnetic fields, we 
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Fig. 3. — Dissipation per unit mass (dotted) and ratio of dissipation to accretion stress 
(solid), versus height. The stress is scaled by the shear rate |f2 so the ratio is dimensionless. 
Both quantities are averaged horizontally and over time from 50 to 200 orbits. 



have verified that MRI modes (k z 3> k r ) grow fastest near the base, where the field is weak, 
while Parker modes (k r ^> k z ) grow fastest high in the atmosphere, where the plasma beta 
is low and magnetic buoya ncy strong. A similar pa ttern appears in ideal-MHD models of 
strongly-magnetized disks ( Uohansen fc Levin! I2008I ). The Parker modes grow faster than 
the orbital frequency in locations where the fields decline steeply with height, including just 
below reversals in the sign of the dominant toroidal field, and along the tops of magnetic flux 
tubes. Among the solutions of the dispersion relation are still faster-growing Parker modes 
with azimuthal wavelengths 2-3 times our domain size. Future experiments with larger boxes 
are indicated. Since the dissipation is concentrated in thin current sheets (figure H]), it dom- 
inates the heating at some places and times. Dissipation is stronger than irradiation in a 
volume fraction 3.2 x 10 -4 , time- averaging from 50 to 200 orbits and considering the layers 
above the disk photosphere (\z\ > 0.18 AU). Temperatures in the current sheets are up to 
50% greater than in the gas heated only by starlight. Finally, mass leaves the top and bot 



torn b oundaries in an outflow comparable to those in isothermal calculations by ISuzuki et al. 

fcoioh . 

Most of the energy input between 50 and 200 orbits (97%) comes from stellar irradiation. 
Of the 3% coming from work done by magnetic stresses, 52% is converted directly to heat 
via magnetic dissipation, while 26% is first converted to kinetic energy and then dissipated 
via shocks and viscosity. The remainder crosses the boundaries as magnetic and turbulent 
kinetic energy and the gravitational potential energy of the small amount of escaping gas. 
Overall, energy ultimately leaves the box almost entirely by radiation diffusion. 



Fig. 4. — Left: Typical snapshot of the magnetic field lines, made at 155 orbits. The field 
strength in Gauss is indicated by colors: stronger fields are red, weaker fields blue. Inside the 
yellow isosurfaces the magnetic dissipation rate exceeds 10 -9 erg cm -3 s -1 . Field lines drawn 
in green are those passing through a zone of strong magnetic dissipation at top center. Right: 
close-up of the green field lines, viewed from the star (top) and along the orbit (bottom). 
The dissipation is concentrated in a current sheet of narrow radial extent, separating two 
counter-undulating bundles of field lines. 



The surface where the starlight is absorbed moves up and down irregularly over intervals 
of several orbits (Figure [5]) as fluctuations in the magnetic fields lead to changes in the density 
profile. The fluctuations are associated with episodic reversals in the toroidal field, which 
lift away from the midplane at a rate that increases with height, producing the butterfly- 
wing pattern visible in figure EJ A majority of the volume-integrated accretion heating 
occurs near these reversals. As in pr e vious stratified MRI turbulence c alculations (e.g. 



Miller & Stondl2000l : iHirose et al.ll2006l : Blaes et all 120071 iFlaig et al.ll2010[ ). the density in 



the magnetically-dominated atmosphere fluctuates horizontally by more than a decade, and 
the velocity dispersion exceeds the sound speed, approaching the Alfven speed. 



4. DISCUSSION 

Diagnosing conditions in planet-forming disks requires an understanding of the power 
sources for the observed line and continuum emission. As a step in this direction, we car- 
ried out energy-conserving radiation-MHD calculations of heating and cooling in a local 
shearing box. The starlight that dominates the overall heating is absorbed in an extended, 
magnetically-supported atmosphere. The surface of optical depth unity lies on average 
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Fig. 5. — Accretion stress vs. height and time, together with the surfaces where the starlight 
is absorbed (white lines) and the Rosseland and Planck mean disk photospheres (black dotted 
and solid lines, respectively). The stress is horizontally-averaged and is plotted in units of 
dyn cm -2 . 



1.5 times higher than in the corresponding hydrostatic model. The fraction of the stellar 
luminosity intercepted will be larger in proportion, with consequences for interpreting the 
spectral energy distributions of protostellar disks. The puffy, magnetically-supported atmo- 
spheres may account for the puzzling fact that many young s tars have near-infrared excesse s 



too large to explain with hydrostatic models, as reviewed by lDullemond fc Monnierl ( 120101 ). 
Greater dust mass fractions than our assumed value of 10~ 4 would raise the opacities, push- 
ing the absorbing surface yet higher, while also increasing the recombination cross-section, 
making the weakly-ionized dead zone thicker. 

The accretion flow occurs in the disk atmosphere, not the interior as assumed in clas- 
sical viscous models, and the dissipation occurs still higher (figured]). The viscous model's 
midplane temperature bump vanishes in the MHD case: relocating the dissipation to the 
atmosphere leaves the interior cooler and isothermal (figure |2J). The magnetic fields are gener- 
ated in an optically-thin turbulent layer and rise buoyantly as they dissipate. While the total 
mass column is 1 000 g cm -2 , over half (56%) of the accretion heat is deposited at columns 
within 1 g cm" 2 of the boundaries, and more than one-fifth (23%) within 0.1 g cm -2 . The 
turbulent layer is optically-thin to its own continuum emission independent of the dust deple- 



tion, so long as the grains control both the opacity and the recombination ( )Bai fc Goodman 



20091 ). Because the accretion power heats optically-thin gas, the resulting spectral lines can 
potentially be used to measure the mass flow rate. The dissipation is spatially inhomoge- 
neous, heating gas locally above the irradiation temperature despite the starlight dominating 
the volume-integrated power. Even a small disk patch thus produces emission at a range of 
temperatures. 



Since the disk atmosphere is magnetically-supported, its structure varies in response to 



- 10 - 



changes in the strength of the fields. In particular, the height where the starlight is absorbed 
moves up and down by almost a factor of two over timescales of several orbits, with an 
amplitude around 10% of the distance to the star. The resulting changes in the stellar 
illumination of the disk surface could c ontribute to the widespread low-amplitude infrare d 
variability among young stars with disks ([Morales- Calderon et al.ll2009t iLuhman et al.ll2010l ). 

To make the calculations feasible, we simplified the situation in several respects. Among 
these, we capped the magnetic diffusivity well above the threshold for shutting off the MRI, 
but well below the value computed from the midplane ionization fraction. The weak mid- 
plane magnetic oscillations visible in figure [5] are an artifact of the reduced midplane diffu- 
sivity, and would disappear if the cap were lifted above the threshold where diffusion cancels 
the winding-up of radial fields. The dead zone stresses are unimportant however, as they 
contribute only 2.6% to the total mean mass accretion rate. Also, we treated the stellar irra- 
diation in a horizontally-averaged fashion and kept the starlight incidence angle fixed despite 
the changes in the disk thickness, since the shearing-box approximation gives no informa- 
tion on the tilt of the disk surface. Future work could explore the effects of relaxing these 
simplifications. In addition, our current sheets are artificially broadened to the grid scale by 
the numerical resistivity. The sheets should be narrower, with still greater accretion heating 
per unit volume, and thus greater chances of exceeding the stellar irradiation. The magnetic 
dissipation furthermore heats the gas, which may become hotter than the dust where den- 
sities are low enough so heat transfer by gas-dust collisions is slow. Since we assumed good 
gas-du st thermal coup li ng, ou r peak temperatures a r e low er lim its. Future modeling sho uld 
follow iGlassgold etaD j2004h . iKamp fc Dullemondl (120041 ) and iNomura fc Millarl hoO$ \ in 
allowing distinct dust and gas temperatures. 
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